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QCD spectral functions of hadrons in the pseudo-scalar 
and vector channels are extracted from lattice Monte Carlo 
data of the imaginary time Green's functions. The maximum 
entropy method works well for this purpose, and the reso- 
nance and continuum structures in the spectra are obtained 
in addition to the ground state peaks. 
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Among various dynamical quantities in quantum chro- 
modynamics (QCD), the spectral functions (SPFs) of 
hadrons play a special role in physical observables (See 
e.g. 0J^]). A well-known example is the cross section 
of the e + e _ annihilation into hadrons, which can be ex- 
pressed by SPF in the vector channel. SPF at finite tem- 
perature (T) and/or baryon density is also a key concept 
to understand the medium modification of hadrons ||. 
The enhancement in low-mass dileptons observed in rel- 
ativistic heavy ion collisions at CERN SPS [| is a typ- 
ical example which may indicate a spectral shift in the 
medium Q. 

However, the Monte Carlo simulations of QCD on the 
lattice, which have been successful in measuring static 
observables ||, have difficulties in accessing the dynam- 
ical quantities in the Minkowski space such as SPFs and 
the real time correlation functions. This is because mea- 
surements on the lattice can only be carried out for dis- 
crete points in imaginary time. The analytic continuation 
from the imaginary time to the real time using the finite 
number of lattice data with noise is highly non-trivial 
and is even classified as an ill-posed problem. 

In this paper, we make a first serious attempt to 
extract SPFs of hadrons from lattice QCD data with- 
out making a priori assumptions on the spectral shape. 
For this purpose, we use the maximum entropy method 
(MEM), which has been successfully applied for similar 
problems in quantum Monte Carlo simulations in con- 
densed matter physics, image reconstruction in crystal- 
lography and astrophysics, and so forth @,^|. Due to the 
limitation of space, we present only the results for the 
pseudo-scalar (PS) and vector (V) channels at T = 0. 
The results for other channels will be given in 0]. 

The Euclidean correlation function D(t) of an opera- 
tor 0(t,x) and its spectral decomposition at zero three- 
momentum read 

D(t) = f(O^(T,x)O{0,6))d 3 x= f K{r,uj)A{uj)duj, (1) 



where r > 0, uj is a real frequency, and A(lj) is SPF 
(or sometimes called the image in this paper), which is 
positive semi-definite by definition. The kernel K(t, uj) 
is proportional to the Fourier transform of a free boson 
propagator with mass uj: At T = 0, K(t, uj) — exp(— tlj). 

Monte Carlo simulation provides D(ri) for the discrete 
set of points < n/a < N T , where N T is the temporal 
lattice size and a is the lattice spacing. In the actual 
analysis, we use data points for T m i n < T{ < T max . From 
this data set with noise, we need to reconstruct the con- 
tinuous function A(uj) on the right hand side of (|l]), or 
to make the inverse Laplace transform. This is a typi- 
cal ill-posed problem, where the number of data is much 
smaller than the number of degrees of freedom to be re- 
constructed. This makes the standard likelihood analysis 
and its variants inapplicable jl0| unless strong assump- 
tions on the spectral shape are made. MEM is a method 
to circumvent this difficulty by making a statistical infer- 
ence of the most probable image as well as its reliability 

0- 

The theoretical basis of MEM is the Bayes' theorem in 
probability theory: P[X\Y] = P[Y\X]P[X]/P[Y], where 
P[A|Y] is the conditional probability of X given Y. Let 
D stand for Monte Carlo data with errors for a specific 
channel on the lattice and H summarize all the defini- 
tions and prior knowledge such as A(lo) > 0. The most 
probable image A(uj) for given lattice data is obtained by 
maximizing the conditional probability P[A\DH], which, 
by the Bayes' theorem, is rewritten as 



P[A\DH] cx P[D\AH]P[A\H], 



(2) 



where P[D\AH] (P[A\H]) is called the likelihood func- 
tion (the prior probability). 

For the likelihood function, the central limiting theo- 
rem leads to P[D\AH] = 1 exp(-L) with 



L=\ YX D fc) ~ D A (r l ))C- ] 1 (D(r J ) - D A (r J )), 



(3) 



1,3 



where i and j run from r min /a through T max /a. Z^ is 
a normalization factor given by Zl = (27r) JV / 2 v / det C 
with N = r max /a — T m in /a +1- D(ri) is the lattice 
data averaged over gauge configurations and D A (ri) is 
the correlation function defined by the right hand side of 
(0). C is an N x N covariance matrix defined by dj = 

[N conf (N conf - I)]' 1 YlzV{D m {n) - DWipnfa) - 
D{tj)): Here N con f is the total number of gauge con- 
figurations and D m (ri) is the data for the m-th gauge 
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configuration. The lattice data have generally strong cor- 
relations among different r's, and it is essential to take 
into account the off-diagonal components of C. 

It can be generally shown on an axiomatic basis |TI]] 
that, for positive distributions such as SPF, the prior 
probability can be written with parameters a and m 
as P[A\Ham] = Zg 1 exp(aS'). Here S is the Shannon- 
Jaynes entropy, 



S 



A{ui) - m(uj) - A(w) log 



duj. (4) 



Zs is a normalization factor: Zg = j e aS VA. a is a real 
and positive parameter and m(u>) is a real function called 
the default model. 

In this paper, we adopt a state-of-art MEM j7j, where 
the output image A out is given by a weighted average 
over A and a: 

A out (uj) = J A(lu) P [A | DHam] P [a \DHm] VA da 

A a {uo) P[a\DHm] da. (5) 



Here A a (oj) is obtained by maximizing Q = aS — L for 
a given a, and we assume that P[A\DHam] is sharply 
peaked around A a (uj). a dictates the relative weight of 
the entropy S (which tends to fit A to the default model 
m) and the likelihood function L (which tends to fit A 
to the lattice data). Note, however, that a appears only 
in the intermediate step and is integrated out in the fi- 
nal result. We found that the weight factor P[a\DHm], 
which is calculable using Q [j7j, is highly peaked around 
its maximum a = a in our lattice data. One can also 
study the stability of the A out (uj) against a reasonable 
variation of m(u>). 

The non-trivial part of the MEM analysis is to find 
the global maximum of Q in the functional space of 
A{uS), which has typically 750 degrees of freedom in our 
case. We have utilized the singular value decomposition 
method to define the search direction in this functional 
space. The method works successfully to find the global 
maximum within reasonable iteration steps. The techni- 
cal details will be given in H . 

To check the feasibility of the MEM procedure and to 
see the dependence of the MEM image on the quality of 
the data, we made the following test using mock data, 
(i) We start with an input image Ai n (u>) = uj 2 pi n {us) 
in the p-meson channel which simulates the experimen- 
tal e + e~ cross section. Then we calculate Di n (r) from 
A in (uj) using eq.(|lj). (ii) By taking An fa) at N dis- 
crete points and adding a Gaussian noise, we create a 
mock data Anocfcfa)- The variance of the noise trfa) 
is given by ofa) = b x Anfa) x T i/ a with a parameter 
b, which controls the noise level 112). (iii) We construct 
the output image A out (uj) = Lu 2 p out (uj) using MEM with 
D moc k(T m in < n < T max ) and compare the result with 



Ai n (uj). In this test, we have assumed that C is diagonal 
for simplicity. 

In Fig.l, we show pi n (uS), and p ut{w) for two sets of 
parameters, (I) and (II). As for m, we choose a form 
m(u>) — mow 2 with mo = 0.027, which is motivated 
by the asymptotic behavior of A in perturbative QCD, 
A(oj > lGeV) = (1/4tt 2 )(1 + a s /ir)uj 2 . The final result 
is, however, insensitive to the variation of mo even by fac- 
tor 5 or 1/5. The calculation of A out (uj) has been done 
by discretizing the w-space with an equal separation of 
10 MeV between adjacent points. This number is chosen 
for the reason we shall discuss below. The comparison 
of the dashed line (set (I)) and the dash-dotted line (set 
(II)) shows that increasing r rnax and reducing the noise 
level b lead to better SPFs closer to the input SPF. 

We have also checked that MEM can nicely reproduce 
other forms of the mock SPFs. In particular, MEM works 
very well to reproduce not only the broad structure but 
also the sharp peaks close to the delta-function as far as 
the noise level is sufficiently small. 

We have then applied MEM to actual lattice data. 
For this purpose, quenched lattice QCD simulations have 
been done with the plaquette gluon action and the Wil- 
son quark action by the open MILC code with minor 
modifications f0|| . The lattice size is 20 3 x 24 with 
(3 = 6.0, which corresponds to a — 0.0847 fm (a -1 = 2.33 
GeV), k c = 0.1571 [0, and the spatial size of the lattice 
L s a = 1.69 fm. Gauge configurations are generated by 
the heat-bath and over-relaxation algorithms with a ra- 
tio 1 : 4. Each configuration is separated by 1000 sweeps. 
Hopping parameters are chosen to be k = 0.153, 0.1545, 
and 0.1557 with N con f = 161 for each k. For the quark 
propagator, the Dirichlet (periodic) boundary condition 
is employed for the temporal (spatial) direction. To cal- 
culate the two-point correlation functions, we adopt a 
point-source at x = and a point-sink averaged over 
the spatial lattice-points to extract physical states with 
vanishing three-momentum. For the PS and V channels, 
the operators dj^u and dj^u (p — 1, 2, 3) are chosen, 
respectively. We use data at 1 < Ti/a < 12 to remove 
the noise at the Dirichlet boundary. To avoid the known 
pathological behavior of the eigenvalues of C Q, we take 

N con f > N. 

We define SPFs for the PS and V channels as 



A(co) 



uj 2 p r 



•M, 



(6) 



so that p PS v (lo — > large) approaches a finite constant as 
predicted by perturbative QCD. For the MEM analysis, 
we need to discretize the w-integration in (^). Since Auj 
(the mesh size) « 1/ r max should be satisfied to suppress 
the discretization error, we take Acu = 10 MeV. Lu m ax (the 
upper limit for the u> integration) should be compara- 
ble to the maximum available momentum on the lattice: 
u m ax ~ tt/o ~ 7.3 GeV. We have checked that larger 
values of uj m ax do not change the result of A{oj) sub- 
stantially, while smaller values of io max distort the high 
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energy end of the spectrum. The dimension of the image 
to be reconstructed is = uj max /Auj ~ 750, which is 
in fact much larger than the maximum number of Monte 
Carlo data N = 25. In the MEM analysis presented in 
this paper, the continuum kernel K — exp(— Ttj) is used. 
We have also carried out analysis based on the lattice 
kernel K lat . A comparison of the two cases will be given 
in j|] in detail. 

In Fig. 2 (a) and (b), we show the reconstructed images 
for each k. In these figures, we have used to = mow 2 with 
too = 2.0(0.86) for PS (V) channel motivated by the per- 
turbative estimate of Too (see eq.(|7]) and the text below). 
We have checked that the result is not sensitive, within 
the statistical significance of the image, to the variation of 
Too by factor 5 or 1/5. The obtained images have a com- 
mon structure: the low-energy peaks corresponding to 7r 
and p, and the broad structure in the high-energy region. 
From the position of the pion peaks in Fig. 2(a), we ex- 
tract k c = 0.1570(3), which is consistent with 0.1571 JIJ] 
determined from the asymptotic behavior of D(t). The 
mass of the p-meson in the chiral limit extracted from 
the peaks in Fig. 2(b) reads m p a — 0.348(15). This is 
also consistent with m p a = 0.331(22) @ determined by 
the asymptotic behavior. Although our maximum value 
of the fitting range r max /a = 12 marginally covers the 
asymptotic limit in r, we can extract reasonable masses 
for 7r and p. The width of 7r and p in Fig. 2 is an artifact 
due to the statistical errors of the lattice data. In fact, 
in the quenched approximation, there is no room for the 
p-meson to decay into two pions. 

As for the second peaks in the PS and V channels, 
the error analysis discussed in Fig. 3 shows that their 
spectral "shape" does not have much statistical sig- 
nificance, although the existence of the non-vanishing 
spectral strength is significant. Under this reservation, 
we fit the position of the second peaks and made lin- 
ear extrapolation to the chiral limit with the results, 
m 2nd /m p = 1.88(8) (2.44(11)) for the PS (V) channel. 
These numbers should be compared with the experimen- 
tal values: m K naao)/ m p = 1.68, and TO p ( 1 45 )/m p = 1.90 
or ™ p (i700)/ m p = 2 - 2 0. 

One should remark here that, in the standard two- 
mass fit of D(t), the mass of the second resonance is 
highly sensitive to the lower limit of the fitting range, 
e.g., m 2nd /m p = 2.21(27) (1.58(26)) for r min /a = 8(9) 
in the V channel with (3 = 6.0 H]. This is because the 
contamination from the short distance contributions from 
t < T m i n is not under control in such an approach. On 
the other hand, MEM does not suffer from this difficulty 
and can utilize the full information down to r m j n /a = 1. 
Therefore, MEM opens a possibility of systematic study 
of higher resonances with lattice QCD data. 

As for the third bumps in Fig. 2, the spectral "shape" is 
statistically not significant as is discussed in Fig. 3, and 
they should rather be considered a part of the pertur- 
bative continuum instead of a single resonance. Fig. 2 



also shows that SPF decreases substantially above 6 GeV; 
MEM automatically detects the existence of the momen- 
tum cutoff on the lattice ~ n/a. It is expected that MEM 
with the data on finer lattices leads to larger ultraviolet 
cut-offs in the spectra. The height of the asymptotic form 
of the spectrum at high energy is estimated as 

p v (^6GeV) = ^(l + ^) y^) 2 -0.86. (7) 

The first two factors are the qq continuum expected 
from perturbative QCD. The third factor contains the 
non-perturbative renormalization constant for the lattice 
composite operator. We adopt Z v = 0.57 determined 
from the two-point functions at j3 = 6.0 [[15j together 
with a s = 0.21 and k = 0.1557. Our estimate in eq.(Q) 
is consistent with the high energy part of the spectrum 
in Fig. 2(b) after averaging over lo. We made a similar 
estimate for the PS channel using Z ps — 0.49 |l^] and 
obtained p PS (w — 6GeV) ~ 2.0. This is also consistent 
with Fig. 2(a). We note here that an independent anal- 
ysis of the imaginary time correlation functions [p| also 
shows that the lattice data at short distance is dominated 
by the perturbative continuum. 

Within the MEM analysis, one can study the statis- 
tical significance of the reconstructed image by the fol- 
lowing procedure 0. Assuming that P[A\DHam] has a 
Gaussian distribution around the most probable image 
A, we estimate the error by the covariance of the image, 
— ({SaSaQ)^ 1 ) a=A' where 5a is a functional derivative 
and (•) is an average over a given energy interval. The fi- 
nal error for A out is obtained by averaging the covariance 
over a with a weight factor P[a\DHm\. Shown in Fig. 3 
is the MEM image in the V channel for k = 0.1557 with 
errors obtained in the above procedure. The height of 
each horizontal bar is (pout(^>)) in each u) interval. The 
vertical bar indicates the error of (p ut(u))- The small 
error for the lowest peak in Fig. 3 supports our identifi- 
cation of the peak with p. Although the existence of the 
non-vanishing spectral strength of the 2nd peak and 3rd 
bump is statistically significant, their spectral "shape" is 
either marginal or insignificant. Lattice data with better 
quality are called for to obtain better SPFs. 

In summary, we have made a first serious attempt to 
reconstruct SPFs of hadrons from lattice QCD data. We 
have used MEM, which allows us to study SPFs without 
making a priori assumption on the spectral shape. The 
method works well for the mock data. Even for the lattice 
data, the method produces resonance and continuum-like 
structures in addition to the ground state peaks. The sta- 
tistical significance of the image has been also analyzed. 
Better data with finer and larger lattice will produce bet- 
ter images with smaller errors, and our study should be 
considered a first attempt towards this goal. We have 
not made the chiral extrapolation of SPFs in this paper, 
since we have found that neither the direct extrapolation 
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of the MEM image A out (uj) nor the extrapolation of D(t) 
and C works in a straightforward manner. We leave this 
as an open problem. From the physics point of view, the 
spectral change at finite temperature is by far the impor- 
tant problem. This is currently under investigation. 
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(No. 10874042) of the Japanese Ministry of Education, 
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FIG. 1. The solid line IS pin (lu). The dashed line and 
dash-dotted line are p ut(u) obtained with parameter set (I) 
a = 0.0847 fm, 1 < r/a < 12, b = 0.001 and set (II) 
a = 0.0847 fm, 1 < r/a < 36, b = 0.0001, respectively. 
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FIG. 2. Reconstructed image p ut(u) for the PS (a) and 
V (b) channels. The solid, dashed, and dash-dotted lines are 
for k = 0.1557, 0.1545, and 0.153, respectively. For the PS 
(V) channel, mo is taken to be 2.0 (0.86). u max is 7.5 GeV in 
this figure and Fig. 3. 




FIG. 3. pout(oj) in the V channel for k = 0.1557 with error 
attached. 
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